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Abstract 

Bridging algorithms are global Monte Carlo moves which allow for an efficient sam- 
pling of single polymer chains. In this manuscript we discuss the adaptation of three 
bridging algorithms from lattice to continuum models, and give details on the correc- 
tions to the acceptance rules which are required to fulfill detailed balance. For the first 
time we are able to compare the efficiency of the moves by analyzing the occurrence 
of knots in globular states. For a flexible homopolymer chain of length = 1000, 
independent configurations can be generated up to two orders of magnitude faster than 
with slithering snake moves. 

Key words: off-lattice Monte Carlo simulations, globular polymers, knots 



1. Introduction 

The estimation of thermodynamic quantities for globular phases of single chain 
polymers remains a challenging task Uj 2, 3] because an effective sampling is hindered 
by self-entanglements, knots flSS and high density. Globular polymers display 
an interesting phase behavior 1 9lll0ll a nd may serve as model systems to study aspects 
of biological macromolecules 1 llll2Lri3. 141 . Typically, approaches like the slithering 
snake algorithm ([l5 , 3, 12 , 18 1 or various types of chain growth algorithms like PERM 
ifigll are applied to sample such configurations. 

A new class of Monte Carlo moves which takes advantage of destroying and re- 
connecting bonds was first suggested independently by Mansfield li20ll and Olaj and 
Lantschbauer fl^ in 1982. Even though both papers focus on sampling (polydisperse) 
lattice melts, Ref. [20] also includes the so called backbite move which preserves 
monodispersity and can be applied to single chains as well. Two additional bridg- 
ing moves for single polymer chains on the lattice were subsequently introduced by 
Deutsch in 1997 |22.] and have since been combined with schemes to sample free en- 



ergy landscapes 11 1 3L 1911 Uke Wang-Landau sampling II23I1 . Nowadays, a rather evolved 
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set of bridging methods exists for off-lattice atomistic polymer melts |24, 251 2^ 27, 
2811290, and one particular move was also ported to single chain simulations ISOll . Apart 



from the obvious application of sampling static melt quantities, these moves are also 
applied to generate well-equilibrated starting configurations for subsequent Molecular 
Dynamics runs ll3lll32|] . 

Our paper focuses on the efficient generation of globular off-lattice polymers. To 
this extent we adapt the three single chain bridging moves of Mansfield and Deutsch 
1 2211 to continuum models. Technical aspects of our implementation which further im- 
prove efficiency are also discussed. In the second part we test the performance of the 
moves and compare it to a standard implementation of the slithering snake algorithm. 
Distinguishing independent configurations in globular states is a formidable task on its 
own as commonly used observables are typically unable to yield information about the 
topology of the chain. To this end, we determine whether or not a globular state is knot- 
ted and compare the times which are required to generate "topologically independent" 
configurations. 



2. Description of the bridging moves and detailed balance 

2.1. Backbite move 

On a lattice, the number of neighbors of a particular monomer is always well- 
defined. In the continuum, we call two monomers neighbors if the distance d between 
them is within a certain range d„j„ < d < dmax- The basic idea of the backbite move is 
illustrated in figure [T] 

First, we select an end monomer with equal probabiUty. Neighbors of this monomer 
are identified and counted (rioid^nd), and one of these neighbors is chosen at random. 
(Fig. [TJ)). A new bond is created between the neighbor and the old end monomer. 
At the same time, the bond between the chosen neighbor and its successor (along the 
direction from the old end to the neighbor) is severed, turning the neighbor into the 
new terminus (Fig.[T];). Note that the selection probability for the reverse Monte Carlo 
move Oyx (from y to x) may be different from a^,.. Therefore, we also need to count the 
number of neighbors of the new end n„e„^„d (Fig.[I]d): 

axy^ and a,,, = - ■ , (1) 

Moves are accepted with a modified Metropolis criterion: 



r < min 



1 , exp (-/3AV) 

^newjsnd 



(2) 



where re (0, 1) is a random number between and 1 . The prefactor ensures that 

new .end 

the asymmetric proposal probability is corrected. AV denotes the energy difference 
between the new state y and the old state x. 

The backbite move was first suggested by Mansfield for a lattice model flB]. Note 
that for a lattice model, the number of neighbors need not be counted necessarily. One 
can simply choose a potential site and reject the move if there is no monomer present. 
Successful implementations in the continuum were reported recently in Refs. iioilsl. 
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2.2. Internal long range move of type II 

In Ref. 1 22] Deutsch proposes two additional moves for single chains in globular states, 
which were originally used to search for ground-states in the HP-model |34]. Both 
Monte Carlo moves cut and "rewire" the polymer internally, and have to our knowledge 
not been ported to a continuum model so far First, we recall Deutsch's definition of 
parallel and anti-parallel bonds; Two opposing strands along the chain made up of 
monomers /, ; + 1 and /', k are named parallel if k-j = -t- 1 and anti-parallel ifk- j - - 1 
(see figs. |2] Ell. In the continuum monomers which are about to be rewired need to be 
located within d„n„ < d < d,„ax. In this notation new bonds are always created between 
i and j and / + 1 and k, which need to be neighbors. When bonds are anti-parallel 
(and not parallel as suggested in the original publication), we attempt a move of type I. 
Likewise, if bonds are parallel we try a type II move. For pedagogical reasons we will 
first discuss moves of type II (Fig.|2]l because they are easier to implement. 

First we select one bond b\ made up of monomers (/, / -H 1) at random (with proba- 
bility l/iN - 1)) and choose one of n/,i neighboring parallel bonds ^2- If no such bonds 
exist, the move is rejected. Bonds bi between / and / + 1 and between b2 between j and 
k are cut and reconnected to bj, between / and j and Z74 between / + 1 and k. This leads 
to following selection probabilities from state x to y: 

1/1 1\ 

+ — ■ (3) 



N -l\nbi ni,. 

Note that the same configuration could have been chosen by first selecting bi instead 
of bi, hence we get the additional summand. For the reverse move, we obtain 

„„ = ^/±.±], ,4, 



- 1 \ Hfc, nt 



A move is accepted if 



r < mm 



nb^nb,(nb, + nbj 

1 , exp {-/3AV) 

nb,nb^{nbt + nh^) 



(5) 



As before, AV denotes the energy difference between the new state y and the old state 
X. Note that the selection probabilities are symmetric in the original publication of 
Deutsch |22]. One simply checks for neighboring bonds on the lattice and rejects the 
move if no adequate bond is present. 



2.3. Internal long range move of type I 

Figure Oprovides details on the long range move of type I. First we select one bond 
at random (with probability 1/{N-1)) and choose one of the neighboring anti- 
parallel bonds b2 (with probability 1/nii). Bonds bi between / and / + 1 and between 
b2 between j and k are cut and reconnected to b^ between / and j and ^74 between / + 1 
and k such that the chain is split into a linear and a circular part. From the circular 
part which consists of bonds we choose a random bond (^5 in fig.[3]with probability 
1 /N^) and check if neighboring bonds belong to the linear part. One of these bonds 
is chosen at random {b(y) and the chain is "rewired" to form the anti-parallel bonds bj 
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and bs- Sometimes it is possible to reconnect bonds ^5 and bd in different ways (as 
indicated by dotted lines in fig.|2]l. In this case each potential connection contributes 
to n^^ . If there is no such bond b(y, the move is rejected. Note that in all steps, the 
number of neighboring bonds need to be determined. This leads to following selection 
probabilities from state x to y: 



1 



1 1 



- 1 n, n' 



1 1 
h — 

"i, tib, 



(6) 



Note that the same configuration could have been chosen by first selecting b2 instead 
of b]_, hence we obtain the additional summand. For the reverse move, we get 



1 



1 1 



A^- 1 



1 1 



This leads to the following Metropolis criterion: 



r < mm 



1, 



nbjtihyi^^inb, + rib,) 



exp (-ySAV) 



(7) 



(8) 



As before, AV denotes the energy difference between the new state y and the old state 

X. 



2.4. Implementation 

In our program monomers of the chain are stored in sequence in a double-linked 
list. This allows for an efficient detection of parallel and anti -parallel bonds for internal 
moves of type II and I. Even more importantly, we set up neighbor tables in which the 
neighbors of each monomer are stored. Unfortunately, this neighbor table needs to be 
updated each time local moves, which are essential for the ergodicity of the algorithms, 
are attempted. In order to minimize these efforts, we always perform several local 
displacements for all particles before we return to a sequence of bridging moves. 

As indicated in previous sections, the internal move of type I is by far the most 
complicated to implement because we need to take care of the chain separation in both 
linear and circular parts. The easiest move is implement is the backbite move followed 
by the internal bridging move of type II. 



3. Model and performance analysis 



In the following we compare the performance of the moveset with an implemen- 
tation of the slithering-snake algorithm UtII : We choose with probability 0.5 one end 
of the chain and attempt to attach the end monomer to a random position at the op- 
posite terminus without changing the bond length. Our model system consists of a 
simple flexible Lennard-Jones + FENE homopolymer Issl 0, [Si with cut and shifted 
Lennard-Jones interactions between all monomers 



4e 



, otherwise 



(9) 
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and FENE interactions between adjacent beads: 



ypENE - -33.75eln 



1 



(l.5cr) 



(10) 



Note that for a homopolymer as described in Eqs. |9]and[T0l the energy difference in 
the Metropolis criterion of the bridging moves only consists of changes in the bond 
energy as the position of the particles are only altered by local moves. All simulations 
took place aiT - 1.66 elkg, which corresponds to 0.57© (for a chain of infinite length 



1 3711.) Therefore, for = 200 to N - 1000 only globular states are observed. 

In an attempt to gauge the efficiency of the three bridging moves, we determined 
whether or not a particular globule is unknotted by calculating the Alexander poly- 
nomial ^2^- Note that mathematically knots are only well-defined in closed curves. 
However, if we connect the ends of an open chain in a well-defined manner we can still 



obtain corresponding information about the "knottedness" of the globule. Ref. II38 I 
provides details on this topic as well as on our implementation of the Alexander poly- 
nomial. The particular closure used in this work is described in Ref. | fl4| ]. To quantify 
the number of independent configurations generated by the bridging algorithms we 
measured the correlation times between unknotted configurations L39il : 



where 



1 0, if the ith configuration contains a knot 
Ui^\ ^ (12) 

1, if the ith configuration is unknotted. 



An "independent configuration" is defined by the decrease of the auto-correlation func- 
tion to e ' . Note that our approach is particularly well-suited for the task of comparing 
bridging moves as it focuses on changes in topology and self-entanglements. In our 
current implementation neither bridging moves nor slithering snake moves are ergodic 
because they do not alter bond lengths. Therefore, it becomes necessary to combine 
the algorithms with local Monte Carlo moves: In our case, each simulation run spent 
26 % of the time with local displacements. 

Fig. m shows the percentage of unknotted configurations for chain lengths ranging 
from = 200 up to = 1000: For = 200 roughly 90 % of all configurations 
are unknotted, for A^ = 1000 only 20 % are unknotted [4]. Figure [5^ shows the num- 
ber of generated configurations per minute divided by our estimate for the geometric 
correlation time as a function of chain length. All configurations were created on a 
single core of a Core2Quad Q6700 processor running at 3.33 Ghz. Figure [Sj? shows 
the same data normalized by the results for the slithering snake moves. Note that the 
time-limiting step in this analysis is not the generation of globular states but the sub- 
sequent knot analysis which took place on a supercomputer. For chains larger than 
A^ = 400 all bridging moves become more efficient than slithering snake moves. The 
internal bridging move of type II is the most efficient followed by the backbite move 
and the internal move of type I. If we consider that the latter is by far the most difficult 
to code among the three, we do not recommend its implementation. Interestingly, the 
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time to generate an "independent configuration" appears to be almost independent of 
chain length for type II and backbite moves, whereas it increases rapidly for slithering 
snake moves as expected. Finally, it is worth noting that combinations of moves may 
even yield better results. A combination of internal II and backbite moves is about two 
orders of magnitude faster for - 1000 than slithering snake moves. 

4. Conclusion 

The aim of this paper is threefold. First, we have implemented three efficient lattice 
bridging moves for single flexible globular homopolymer chains in a continuum model, 
and discussed difficulties arising from the implementation. In an attempt to measure 
performance we determined correlation times between unknotted globular states. To 
our knowledge this is the first time in which topology is considered to gauge the effi- 
ciency of global Monte Carlo moves. From this analysis, we conclude that all bridg- 
ing moves become more efficient than slithering snake moves for chains larger than 

= 400 monomers. Among the three moves, the internal bridging move of type 11 
is the most successful. The most complicated move (internal I) performs worst and 
we do not recommend its implementation. For = 1000 a combination of backbit- 
ing and internal bridging moves of type II is two orders of magnitude faster than our 
implementation of the slithering snake move. 

The exact speed-up factors should, however, be taken with a grain of salt and rather 
serve as a guideline. Long-range correlations, which are only relaxed by local moves, 
were e.g. not considered, and results also depend a lot on the model under consid- 
eration. If chain stiffness is included 1,40.1 , efficiency will drop, and results will also 
deteriorate if the implementation does not include look-up tables. Nonetheless, we be- 
lieve to have demonstrated that bridging algorithms are indeed the method of choice 
for the simulation of long globular polymers in the continuum. 
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Figure 1: Off-lattice version of the backbite move (2oll . 

a) An end monomer is chosen with equal probability. 

b) Neighbors (within distance d„„„ < d < d„,ax) of this monomer are identified and counted. One of these 
neighbors is selected at random. 

c) A new bond is created between the neighbor and the old end monomer The bond between the chosen 
neighbor and its successor (along the direction from the old end to the neighbor) is cut. 

d) The number of neighbors of the new end needs to be counted to fulfill detailed balance. 

e) The chosen neighbor becomes the new end and the new configuration is accepted with a modified Metropo- 
lis criterion. 
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Figure 2: Internal long range move of type n (parallel bonds: k-j=+l) For details see main text. 




Figure 3: Internal long range move of type I (anti-parallel bonds: k-j=-l): The chain is spUt into a linear and 
a circular part and reconnected. For details see main text. 
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Figure 4: Percentage of unknotted configurations for single bead-spring polymer chains of size N=200 to 
N=1000 in the globular phase (T = 1.66 e/ks) 
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Figure 5: Performance analysis for single bead-spring polymer chains of size N=200 to N=1000 in the 
globular phase (T = 1.66 e/kg) 

a. Number of "independent configurations" per minute as defined in the main text. 

b. Speed-up factor of bridging moves in comparison with slithering snake moves. For this particular model 
all rebridging moves become more efficient than sUthering snake moves for chains larger than 400 monomers. 
For clarity we have used logarithmic scale on both axis. 



